23 Juni 2017

Die Pakete ggplot und ggmap

Das Paket ggplot2

  • Entwickelt von Hadley Wickham
  • Viele Informationen unter:
  • http://ggplot2.org/
  • Den Graphiken liegt eine eigene Grammitik zu Grunde

Das Paket ggplot2 installieren und laden

Der diamonds Datensatz

head(diamonds)
carat cut color clarity depth table price x y z
0.23 Ideal E SI2 61.5 55 326 3.95 3.98 2.43
0.21 Premium E SI1 59.8 61 326 3.89 3.84 2.31
0.23 Good E VS1 56.9 65 327 4.05 4.07 2.31
0.29 Premium I VS2 62.4 58 334 4.20 4.23 2.63
0.31 Good J SI2 63.3 58 335 4.34 4.35 2.75
0.24 Very Good J VVS2 62.8 57 336 3.94 3.96 2.48

Wie nutzt man qplot

  • qplot wird für schnelle Graphiken verwendet (quick plots)
  • bei ggplot kann man alles bis ins Detail kontrollieren
# histogram
qplot(depth, data=diamonds)

Ein Balkendiagramm

qplot(cut, depth, data=diamonds)

Ein weiteres Balkendiagramm

qplot(factor(cyl), data=mtcars,geom="bar")

Boxplot

qplot(data=diamonds,x=cut,y=depth,geom="boxplot")

Scatterplot

# scatterplot
qplot(carat, depth, data=diamonds)

Farbe hinzu:

qplot(carat, depth, data=diamonds,color=cut)

Trendlinie hinzufügen

myGG<-qplot(data=diamonds,x=carat,y=depth,color=carat) 
myGG + stat_smooth(method="lm")

Graphik drehen

qplot(factor(cyl), data=mtcars, geom="bar") + 
coord_flip()

Wie nutzt man ggplot

  • die aestetics:
ggplot(diamonds, aes(clarity, fill=cut)) + geom_bar()

Farben selber wählen

Es wird das Paket RColorBrewer verwendet um die Farbpalette zu ändern

install.packages("RColorBrewer")
library(RColorBrewer)
myColors <- brewer.pal(5,"Accent")
names(myColors) <- levels(diamonds$cut)
colScale <- scale_colour_manual(name = "cut",
                                values = myColors)

http://stackoverflow.com/questions/6919025/

Eine Graphik mit den gewählten Farben

p <- ggplot(diamonds,aes(carat, depth,colour = cut)) + 
  geom_point()
p + colScale

Speichern mit ggsave

ggsave("Graphik.jpg")

Links

Verschiedene Kartentypen

Straßenkarten

  • Straßenkarte werden sehr häufig verwendet.
  • Diese Karten zeigen Haupt- und Nebenstraßen (abhängig vom Detail)
  • oft sind auch weitere Informationen enthalten. Wie beispielsweise Flughäfen, Städte, Campingplätze oder andere Orte von Interesse
  • Beispiel einer Straßenkarte für Mannheim.

Installieren des Paketes

  • Zur Erstellung der Karten brauchen wir das Paket ggmap:
devtools::install_github("dkahle/ggmap")
devtools::install_github("hadley/ggplot2")
install.packages("ggmap")

Paket ggmap - Hallo Welt

  • Um das Paket zu laden verwenden wir den Befehl library
library(ggmap)

Und schon kann die erste Karte erstellt werden:

qmap("Mannheim")

Karte für eine Sehenswürdigkeit

BBT <- qmap("Berlin Brandenburger Tor")
BBT

Karte für einen ganzen Staat

qmap("Germany")

  • Wir brauchen ein anderes zoom level

Ein anderes zoom level

  • level 3 - Kontinent
  • level 10 - Stadt
  • level 21 - Gebäude
qmap("Germany", zoom = 6)

Hilfe bekommen wir mit dem Fragezeichen

?qmap

Verschiedene Abschnitte in der Hilfe:

  • Description
  • Usage
  • Arguments
  • Value
  • Author(s)
  • See Also
  • Examples

Die Beispiele in der Hilfe

Ausschnitt aus der Hilfe Seite zum Befehl qmap:

qmap Example

Das Beispiel kann man direkt in die Konsole kopieren:

# qmap("baylor university")
qmap("baylor university", zoom = 14)

# und so weiter

Ein anderes zoom level

qmap("Mannheim", zoom = 12)

Näher rankommen

qmap('Mannheim', zoom = 13)

Ganz nah dran

qmap('Mannheim', zoom = 20)

ggmap - maptype satellite

qmap('Mannheim', zoom = 14, maptype="satellite")

ggmap - maptype satellite zoom 20

qmap('Mannheim', zoom = 20, maptype="hybrid")

ggmap - maptype hybrid

qmap("Mannheim", zoom = 14, maptype="hybrid")

Terrain/physical maps

  • Aus Physischen Karten kann man Informationen über Berge, Flüsse und Seen ablesen.

  • Farben werden oft genutzt um Höhenunterschiede zu visualisieren

ggmap - terrain map

qmap('Schriesheim', zoom = 14,maptype="terrain")

Abstrahierte Karten)

New York

  • Abstraktion wird genutzt um nur essentielle Informationen zu zeigen.
  • Bsp. U-Bahn Karten - wichtig sind Richtungen und wenig Infos zur Orientierung
  • Nun kommen Karten, die sich als Hintergrund eignen.

ggmap - maptype watercolor

qmap('Mannheim', zoom = 14,maptype="watercolor",source="stamen")

ggmap - source stamen

qmap('Mannheim', zoom = 14,
 maptype="toner",source="stamen")

ggmap - maptype toner-lite

qmap('Mannheim', zoom = 14,
 maptype="toner-lite",source="stamen")

ggmap - maptype toner-hybrid

qmap('Mannheim', zoom = 14,
 maptype="toner-hybrid",source="stamen")

ggmap - maptype terrain-lines

qmap('Mannheim', zoom = 14,
 maptype="terrain-lines",source="stamen")

Graphiken speichern

RstudioExport

ggmap - ein Objekt erzeugen

  • <- ist der Zuweisungspfeil um ein Objekt zu erzeugen
  • Dieses Vorgehen macht bspw. Sinn, wenn mehrere Karten nebeneinander gebraucht werden.
MA_map <- qmap('Mannheim', 
               zoom = 14,
               maptype="toner",
               source="stamen")

Geokodierung

Geocoding (…) uses a description of a location, most typically a postal address or place name, to find geographic coordinates from spatial reference data …

Wikipedia - Geocoding

library(ggmap)
geocode("Mannheim",source="google")
lon lat
8.463182 49.48608

Latitude und Longitude

Koordinaten verschiedener Orte in Deutschland

cities lon lat
Hamburg 9.993682 53.55108
Koeln 6.960279 50.93753
Dresden 13.737262 51.05041
Muenchen 11.581981 48.13513

Reverse Geokodierung

Reverse geocoding is the process of back (reverse) coding of a point location (latitude, longitude) to a readable address or place name. This permits the identification of nearby street addresses, places, and/or areal subdivisions such as neighbourhoods, county, state, or country.

Quelle: Wikipedia

revgeocode(c(48,8))
## [1] "Unnamed Road, Somalia"

Die Distanz zwischen zwei Punkten

mapdist("Q1, 4 Mannheim","B2, 1 Mannheim")
##             from             to   m    km     miles seconds  minutes
## 1 Q1, 4 Mannheim B2, 1 Mannheim 749 0.749 0.4654286     215 3.583333
##        hours
## 1 0.05972222
mapdist("Q1, 4 Mannheim","B2, 1 Mannheim",mode="walking")
##             from             to   m    km     miles seconds minutes  hours
## 1 Q1, 4 Mannheim B2, 1 Mannheim 546 0.546 0.3392844     423    7.05 0.1175

Eine andere Distanz bekommen

mapdist("Q1, 4 Mannheim","B2, 1 Mannheim",mode="bicycling")
##             from             to   m    km    miles seconds  minutes
## 1 Q1, 4 Mannheim B2, 1 Mannheim 555 0.555 0.344877     215 3.583333
##        hours
## 1 0.05972222

Geokodierung - verschiedene Punkte von Interesse

POI1 <- geocode("B2, 1 Mannheim",source="google")
POI2 <- geocode("Hbf Mannheim",source="google")
POI3 <- geocode("Mannheim, Friedrichsplatz",source="google")
ListPOI <-rbind(POI1,POI2,POI3)
POI1;POI2;POI3
##        lon      lat
## 1 8.462844 49.48569
##        lon      lat
## 1 8.469879 49.47972
##        lon      lat
## 1 8.475754 49.48304

Punkte in der Karte

MA_map +
geom_point(aes(x = lon, y = lat),
data = ListPOI)

Punkte in der Karte

MA_map +
geom_point(aes(x = lon, y = lat),col="red",
data = ListPOI)

ggmap - verschiedene Farben

ListPOI$color <- c("A","B","C")
MA_map +
geom_point(aes(x = lon, y = lat,col=color),
data = ListPOI)

ggmap - größere Punkte

ListPOI$size <- c(10,20,30)
MA_map +
geom_point(aes(x = lon, y = lat,col=color,size=size),
data = ListPOI)

Eine Route von Google maps bekommen

Eine Karte mit dieser Information zeichnen

qmap("Mannheim Hbf", zoom = 14) +
  geom_path(
    aes(x = lon, y = lat),  colour = "red", size = 1.5,
    data = route_df, lineend = "round"
  )

Wie fügt man Punkte hinzu

Bubbles auf einer Karte

Cheatsheet

Resourcen und Literatur

Take Home Message

Was klar sein sollte:

  • Wie man eine schnelle Karte erzeugt
  • Wie man geokodiert
  • Wie man eine Distanz berechnet

Die lineare Regression

Die lineare Regression

Maindonald - Data Analysis

  • Einführung in R
  • Datenanalyse
  • Statistische Modelle
  • Inferenzkonzepte
  • Regression mit einem Prädiktor
  • Multiple lineare Regression
  • Ausweitung des linearen Modells

Lineare Regression in R - Beispieldatensatz

Das lineare Regressionsmodell in R

Schätzen eines Regressionsmodells:

roller.lm <- lm(depression ~ weight, data = roller)

So bekommt man die Schätzwerte:

summary(roller.lm)
## 
## Call:
## lm(formula = depression ~ weight, data = roller)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.180 -5.580 -1.346  5.920  8.020 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -2.0871     4.7543  -0.439  0.67227   
## weight        2.6667     0.7002   3.808  0.00518 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.735 on 8 degrees of freedom
## Multiple R-squared:  0.6445, Adjusted R-squared:  0.6001 
## F-statistic:  14.5 on 1 and 8 DF,  p-value: 0.005175

Falls das Modell ohne Intercept gesch?tzt werden soll:

lm(depression ~ -1 + weight, data = roller)
## 
## Call:
## lm(formula = depression ~ -1 + weight, data = roller)
## 
## Coefficients:
## weight  
##  2.392

Summary des Modells

summary(roller.lm)
## 
## Call:
## lm(formula = depression ~ weight, data = roller)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.180 -5.580 -1.346  5.920  8.020 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -2.0871     4.7543  -0.439  0.67227   
## weight        2.6667     0.7002   3.808  0.00518 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.735 on 8 degrees of freedom
## Multiple R-squared:  0.6445, Adjusted R-squared:  0.6001 
## F-statistic:  14.5 on 1 and 8 DF,  p-value: 0.005175

R arbeitet mit Objekten

  • roller.lm ist nun ein spezielles Regressions-Objekt
  • Auf dieses Objekt können nun verschiedene Funktionen angewendet werden
predict(roller.lm) # Vorhersage
##         1         2         3         4         5         6         7 
##  2.979669  6.179765  6.713114 10.713233 12.046606 14.180002 14.980026 
##         8         9        10 
## 18.180121 24.046962 30.980502
resid(roller.lm) # Residuen
##          1          2          3          4          5          6 
## -0.9796695 -5.1797646 -1.7131138 -5.7132327  7.9533944  5.8199976 
##          7          8          9         10 
##  8.0199738 -8.1801213  5.9530377 -5.9805017

Residuenplot

  • Sind Annahmen des linearen Regressionsmodells verletzt?
  • Dies ist der Fall, wenn ein Muster abweichend von einer Linie zu erkennen ist.
  • Hier ist der Datensatz sehr klein
plot(roller.lm,1)

Residuenplot

  • Wenn die Residuen normalverteilt sind sollten sie auf einer Linie liegen.
plot(roller.lm,2)

Regressionsdiagnostik mit Basis-R

Ein einfaches Modell

N <- 5
x1 <- rnorm(N)
y <- runif(N)

Die Dichte der beiden Vektoren

par(mfrow=c(1,2))
plot(density(x1))
plot(density(y))

Modellvorhersage machen

mod1 <- lm(y~x1)
pre <- predict(mod1)
y
## [1] 0.08132782 0.25203137 0.02478898 0.67305217 0.51547696
pre
##          1          2          3          4          5 
## 0.41571144 0.18834357 0.03351329 0.45645742 0.45265156

Regressionsdiagnostik mit Basis-R

plot(x1,y)
abline(mod1)
segments(x1, y, x1, pre, col="red")

Beispieldaten Luftqualität

library(datasets)
?airquality

Das visreg-Paket

Ein Modell wird auf dem airquality Datensatz geschätzt

install.packages("visreg")
library(visreg)
fit <- lm(Ozone ~ Solar.R + Wind + Temp, data = airquality)
summary(fit)
## 
## Call:
## lm(formula = Ozone ~ Solar.R + Wind + Temp, data = airquality)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.485 -14.219  -3.551  10.097  95.619 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -64.34208   23.05472  -2.791  0.00623 ** 
## Solar.R       0.05982    0.02319   2.580  0.01124 *  
## Wind         -3.33359    0.65441  -5.094 1.52e-06 ***
## Temp          1.65209    0.25353   6.516 2.42e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.18 on 107 degrees of freedom
##   (42 observations deleted due to missingness)
## Multiple R-squared:  0.6059, Adjusted R-squared:  0.5948 
## F-statistic: 54.83 on 3 and 107 DF,  p-value: < 2.2e-16

Visualisierung

par(mfrow=c(2,2))
visreg(fit)

Und dann mit visreg visualisiert.

  • Zweites Argument - Spezifikation erklärende Variable für Visualisierung
visreg(fit, "Wind", type = "contrast")

Visualisierung mit dem Paket visreg

visreg(fit, "Wind", type = "contrast")

Das visreg-Paket

  • Das Default-Argument für type ist conditional.
visreg(fit, "Wind", type = "conditional")

Regression mit Faktoren

Mit visreg können die Effekte bei Faktoren visualisiert werden.

airquality$Heat <- cut(airquality$Temp, 3, 
    labels=c("Cool", "Mild", "Hot"))
fit.heat <- lm(Ozone ~ Solar.R + Wind + Heat, 
    data = airquality)
summary(fit.heat)
## 
## Call:
## lm(formula = Ozone ~ Solar.R + Wind + Heat, data = airquality)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.473 -12.794  -2.686   8.461 107.035 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 48.27682    8.83072   5.467 3.07e-07 ***
## Solar.R      0.06524    0.02145   3.042  0.00296 ** 
## Wind        -3.49803    0.59042  -5.925 3.94e-08 ***
## HeatMild     9.05367    4.88257   1.854  0.06648 .  
## HeatHot     43.13970    5.98618   7.207 8.79e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.9 on 106 degrees of freedom
##   (42 observations deleted due to missingness)
## Multiple R-squared:  0.6554, Adjusted R-squared:  0.6424 
## F-statistic:  50.4 on 4 and 106 DF,  p-value: < 2.2e-16

Effekte von Faktoren

par(mfrow=c(1,2))
visreg(fit.heat, "Heat", type = "contrast")
visreg(fit.heat, "Heat", type = "conditional")

Das Paket visreg - Interaktionen

airquality$Heat <- cut(airquality$Temp, 3, 
labels=c("Cool", "Mild", "Hot"))
fit <- lm(Ozone ~ Solar.R + Wind * Heat, data = airquality)
summary(fit)
## 
## Call:
## lm(formula = Ozone ~ Solar.R + Wind * Heat, data = airquality)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.472 -11.640  -1.919   7.403 102.428 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.48042   17.38219   0.258 0.797102    
## Solar.R        0.07634    0.02137   3.572 0.000538 ***
## Wind           0.05854    1.34860   0.043 0.965458    
## HeatMild      56.72928   18.53110   3.061 0.002805 ** 
## HeatHot       94.68619   18.61619   5.086 1.63e-06 ***
## Wind:HeatMild -4.11933    1.57108  -2.622 0.010054 *  
## Wind:HeatHot  -4.88125    1.74358  -2.800 0.006101 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.28 on 104 degrees of freedom
##   (42 observations deleted due to missingness)
## Multiple R-squared:  0.6825, Adjusted R-squared:  0.6642 
## F-statistic: 37.26 on 6 and 104 DF,  p-value: < 2.2e-16

Steuern der Graphikausgabe mittels layout

visreg(fit, "Wind", by = "Heat",layout=c(3,1))

Das Paket visreg - Interaktionen overlay

fit<-lm(Ozone~Solar.R+Wind*Heat,data=airquality)
visreg(fit,"Wind",by="Heat",overlay=TRUE,partial=FALSE)

Das Paket visreg - visreg2d

fit2<-lm(Ozone~Solar.R+Wind*Temp,data=airquality)
visreg2d(fit2,"Wind","Temp",plot.type="image")

Das Paket visreg - surface

visreg2d(fit2, "Wind", "Temp", plot.type = "persp")

Linkliste - lineare Regression

Die logistische Regression

Agresti - Categorical Data Analysis (2002)

  • Sehr intiutiv geschriebenes Buch
  • Sehr ausführliches begleitendes Skript von Thompson
  • Das Skript eignet sich um die kategoriale Datenanalyse nachzuvollziehen

Faraway Bücher zu Regression in R

Binäre AVs mit glm

  • Die logistische Regression gehört zur Klasse der generalisierten linearen Modelle (GLM)
  • Die Funktion zur Schätzung eines Modells dieser Klasse in heißt glm()
  • glm() muss 1. ein Formel-Objekt mitgegeben werden und 2. die Klasse (binomial, gaussian, Gamma) samt link-Funktion (logit, probit, cauchit, log, cloglog)

Beispieldaten für die logistische Regression

install.packages("HSAUR")
library("HSAUR")
data("plasma", package = "HSAUR")

Logistische Regression mit R

Logistische Regression mit R

plasma_glm_1 <- glm(ESR ~ fibrinogen, data = plasma, 
                    family = binomial())

Weitere Beispieldaten

install.packages("faraway")
library("faraway")
data(orings)
temp damage
53 5
57 1
58 1

Generalisierte Regression mit R - weitere Funktionen

  • Logistisches Modell mit Probit-Link:
probitmod <- glm(cbind(damage,6-damage) ~ temp, 
    family=binomial(link=probit), orings)
  • Regression mit Zähldaten:
modp <- glm(Species ~ .,family=poisson,gala)
  • Proportional odds logistic regression im Paket MASS:
library("MASS")
house.plr<-polr(Sat~Infl,weights=Freq,data=housing)

Linkliste - logistische Regression

Mehrebenenmodelle

Wie sehen die Daten aus?

  • Beispiel Mehrebenenstruktur der Daten

Andres Gutierrez - Multilevel Modeling of Educational Data using R

  • Lineare Modelle erkennen den Cluster-Effekt aufgrund der Intraklassen Korrelation nicht

Beispiel Mehrebenenmodelle

Untersuchungsgegenstand

  • Es sollen die Kenntnisse (Fähigkeiten) von Grundschülern in Mathematik gemessen werden.
  • Dazu werden in einem Schulbezirk zunächst Schulen ausgewählt und anschließend Klassen.
  • Innerhalb der Klassen soll schließlich jeweils eine Stichprobe gezogen und diese getestet werden.

  • Geht man zunächst von einer zufälligen Auswahl von Klassen aus, dann ist die Level-1-Variation durch die Schüler und die Level-2-Variation durch die Klassen bestimmt.

Fragen hierzu

  • Wie wäre die Auswahl der Schulen zu berücksichtigen?

  • Wie kann zusätzlich eine Unterscheidung nach privaten und staatlichen Schulen in die Modellierung eingebracht werden?

Beispiel in Goldstein (2010), Kapitel 1.2

Evaluierung der Effektivität von Schulen

Mehrebenen-Modelle:

  • Schüler
  • Klassenverbände
  • Schulamtsbezirke oder Bundesländer

Unterscheidung

  • Modelle mit vielen Parametern, die wiederum modelliert werden können
  • Regressionen mit Koeffizienten, die zwischen Gruppen variieren können

Bibliotheken

# Linear Mixed-Effects Models using 'Eigen' and S4
install.packages("lme4")

# Data Visualization for Statistics in Social Science
install.packages("sjPlot")
  • Nötige Pakete werden geladen
library(ggplot2)

# Miscellaneous Functions for "Grid" Graphics
library(gridExtra)
library(lme4)
library(sjPlot)

# A Grammar of Data Manipulation
library(dplyr)

Beispieldaten

mlexdat <- read.csv(
"https://github.com/Japhilko/RSocialScience/
raw/master/data/mlexdat.csv") 
X SES Score ID
1 18.62733 -55.120574 A
2 33.64915 -92.375273 A
3 22.26931 -48.783447 A
4 36.49052 38.099329 A
5 38.21402 339.701754 A
6 11.36669 2.286978 A

Formalistisch

  • Bei der Analyse von Daten mit diesen hierarchischen Strukturen, sollte man immer zunächst ein Null-Modell anpassen
  • Somit kann man die Variation erfassen, die auf die Schulen zurückzuführen ist.

  • Das passende Modell sieht folgendermaßen aus:

Die Gesamtvariation wird in zwei Teile zerlegt:

  • Variation zwischen Schülern (innerhalb der Schulen) und
  • zwischen den Schulen (zwischen den Schulen).

Der R-code für dieses Nullmodell

  • das einfachste Multilevel Modell
  • nach dem vertikalen Strich wird die Gruppen Variable spezifiziert
  • die Default Schätzmethode ist restricted maximum likelihood (REML)
  • Man kann aber auch maximum likelihood Schätzung spezifizieren (ML)
HLM0 <- lmer(Score ~ (1 | ID), data = mlexdat)

Nullmodell Ergebnis

coef(HLM0)
## $ID
##   (Intercept)
## A     45.7893
## B    430.7218
## C   1182.1760
## D   2145.2329
## E   3489.1408
## 
## attr(,"class")
## [1] "coef.mer"
summary(HLM0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Score ~ (1 | ID)
##    Data: mlexdat
## 
## REML criterion at convergence: 7130.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.74559 -0.69317 -0.00757  0.68337  2.96511 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 1931758  1389.9  
##  Residual               87346   295.5  
## Number of obs: 500, groups:  ID, 5
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   1458.6      621.7   2.346

Interpratation des Nullmodells

  • 96 Prozent Variation zwischen den Schulen
  • 4 Prozent Variation innerhalb der Schulen
100 * 87346 / (87346 + 1931757)
## [1] 4.32598
  • Die Schätzung der zufälligen Effekte zeigt, dass die Variation zwischen den Schulen (Intraklassen Korrelation) fast 96 Prozent beträgt

  • Während der Anteil der Variation zwischen den Studierenden nur etwas mehr als 4 Prozent ausmacht.

  • Das Null-Modell behauptet also , dass Leistungsträger zu bestimmten Schulen gehen und Studierende mit geringerem Leistungsniveau nicht diese Schulen besuchen.

  • Mit anderen Worten, die Schule bestimmt das Testergebnis.

Ein weiteres Modell

  • Das Null Modell schließt keine erklärenden Variablen ein.
  • Allerdings könnte der sozioökonomischen Status (SES) der Schüler auch eine Rolle spielen.
  • Die folgenden Ausdrücke geben ein verfeinertes Modell mit zufälligen Achsenabschnitten und Steigung für jede der Schulen.

Rcode für dieses Modell

HLM1 <- lmer(Score ~ SES + (SES | ID), data = mlexdat)
coef(HLM1)
## $ID
##   (Intercept)        SES
## A    36.46401  0.3798185
## B    37.21549  9.7596237
## C    38.10719 20.8897245
## D    38.85566 30.2320132
## E    39.70159 40.7907386
## 
## attr(,"class")
## [1] "coef.mer"
summary(HLM1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Score ~ SES + (SES | ID)
##    Data: mlexdat
## 
## REML criterion at convergence: 6742.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.83274 -0.64740  0.02662  0.69063  2.67309 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  ID       (Intercept)     1.65   1.285      
##           SES           257.09  16.034  1.00
##  Residual             40400.24 200.998      
## Number of obs: 500, groups:  ID, 5
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   38.069     45.863   0.830
## SES           20.410      7.236   2.821
## 
## Correlation of Fixed Effects:
##     (Intr)
## SES -0.119

# 1% - BS variance
# 99% - WS variance
100 * 40400.24 / (40400.24 + 257.09 + 1.65)
## [1] 99.36363
# Percentage of variation explained by SES between schools
1 - ((257.09 + 1.65) / 1931757)
## [1] 0.9998661
# Percentage of variation explained by SES within schools
1 - (40400.24 / 87346)
## [1] 0.5374689
  • die Variable SES erklärt 99 Prozent der Unterschiede zwischen den Schulen
  • diese Variable SES erklärt 53 Prozent der Abweichungen innerhalb der Schulen.

Was heißt das? - Schulsegregation

  • wohlhabende Studenten gehören zu reichen Schulen
  • arme Studenten gehören zu armen Schulen.

  • Die Leistung der wohlhabenden Studenten ist höher als die der armen Studenten.

Ein weiteres Beispiel zur Spezifikation von Multilevel Modellen

  • benötigte Bibliotheken:
library(lme4)
library(mlmRev)

Der Datensatz

data(Exam)
# names(Exam)
school normexam schgend schavg vr intake standLRT sex type student
1 0.2613242 mixed 0.1661752 mid 50% bottom 25% 0.6190592 F Mxd 143
1 0.1340672 mixed 0.1661752 mid 50% mid 50% 0.2058022 F Mxd 145
1 -1.7238820 mixed 0.1661752 mid 50% top 25% -1.3645760 M Mxd 142
1 0.9675862 mixed 0.1661752 mid 50% mid 50% 0.2058022 F Mxd 141
1 0.5443412 mixed 0.1661752 mid 50% mid 50% 0.3711052 F Mxd 138
1 1.7348992 mixed 0.1661752 mid 50% bottom 25% 2.1894372 M Mxd 155

Zufälliger Intercept und fixed predictor auf individeller Ebene

  • Ein Prädiktor wird auf jeder Ebene hinzugefügt
  • Dazu wird die '1' im Nullmodell durch den Prädiktor (hier: standLRT) ersetzen.
  • Es wird immer ein Intercept angenommen
  • Da wir nicht wollen, dass der Effekt des Prädiktors zwischen den Gruppen variiert, bleibt die Spezifikation des zufälligen Teils des Modells mit dem vorherigen Modell identisch.
lmer(normexam ~ standLRT + (1 | school), data=Exam)

Random intercept, Random slope

  • Modell mit zufälligen Intercept auf individueller Ebene und
  • einem Prädiktor, der zwischen Gruppen variieren darf.

  • Mit anderen Worten: die Wirkung der Hausaufgaben auf das Ergebnis der Klausur (Mathe-Test) variiert zwischen den Schulen.

  • Zur Schätzung wird '1' - der Intercept im zufälligen Teil der Modellspezifikation
  • …durch die Variable ersetzt, von der wir den Effekt zwischen den Gruppen variieren wollen.

Varying intercept model

MLexamp.6 <- lmer(extro ~ open + agree + social + (1 | school), data = lmm.data)

Varying slope model

MLexamp.9 <- lmer(extro ~ open + agree + social + (1 + open | school/class), data = lmm.data)

Links